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B. D'AURIA 

Abstract. In this paper we investigate an M/M/oo queue whose pa- 
rameters depend on an external random environment that we assume 
to be a semi-Markovian process with finite state space. For this model 
we show a recursive formula that allows to compute all the factorial 
moments for the number of customers in the system in steady state. 
The used technique is based on the calculation of the raw moments of 
the measure of a bidimensional random set. Finally the case when the 
random environment has only two states is deeper analyzed. We obtain 
an explicit formula to compute the above mentioned factorial moments 
when at least one of the two states has sojourn time exponentially dis- 
tributed. 



1. Introduction 

The M /M / oo queue is one of the simplest model in queueing theory. This 
is due to the joint situation to have a memory-less arrival process and an 
infinite set of servers that allows customers to behave independently from 
each other. This suddenly stops to be true after introducing some correla- 
tion between customers. In this paper we achieve that by introducing an 
independent random environment that modulates the system parameters, 
i.e. the arrival rate and the server speeds. Queues with variable service 
and arrival speeds arise naturally in practice and therefore many classi- 
cal works can be foun d. Most of t h e result s deal with the s i ngle s erver 
queue, see for exam ple Takind (120051 ). Ozawa (120041 ). ISeneuptal (|l990l ) and 



references therein. Neutsl ( 198ll ) analyzed the M/M/l queue as well as 



the M/M/C qu eue in random environment by using the matrix-geometric 
approach while Takine and Senguptal ( 1997 ) looked at the infinite server 
queue when only the arrival process was subject to a Markovian modula- 
tio n. The infinite serv e r que u e in random environment has the n bee n studied 
bv iKeilson and Servil (ll993T ). iBavkal-Gursov and Xiad (|2004h and iD'Auria 
(|2007h m the special case when the random environment is Markovian and 
has only two states. 
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In O'Cinneide and Purdue ( 19861 ) the authors looked at the case when the 



environment is given by a finite state Markov process and for this case they 
showed how to compute the factorial moments for the number of customers 
in the system in steady state. Here we extend their analysis to the case of 
a semi-Mar kovian random environment. 

This extension is interesting as it makes the model more attractive for 
application purposes. Indeed, despite its simplicity, the M/M/oo system is 
often used to model pure delay systems, such as highways, satellite links or 
long communication cables, or to approximate the behavior of multi server 
systems. When these kinds of systems are subject to external influences, 
such as day time changing rates, it is then helpful to look at extended models, 
such as the one proposed in this work, in order to analyze or predict their 
behaviors. 



The methodology we use follows the technique developed in iD'Auria 



( 20071 ). It consists on representing the stationary and isolate M/G/oo sys- 



tem as a Poisson process on Mr and by computing the number of customers 
in the system by measuring a det erministic set according to the p o int pr ocess 
measure (see also Resnick ( 1987 . §3.3) and D'Auria and Resnick ( 20061 )). In 



this context the random environment can be expressed as a random modula- 
tion of the set and in the special case of exponential distributed service times 
its measure can be derived by solving a system of stochastic equations, see 
relation (|4.5p below. We use this relation to compute the factorial moments 
for the number of customers in the system at steady state. 

2. Model description 

To start, we define the random environment {T(u),u € R} as a semi- 
Markov chain with values in the finite state space E = {1, . . . ,K}. We 
assume that the sojourn time in the state denoted as Tk, is an inde- 

pendent positive random variable whose distribution function has Laplace 
transform denoted by Tfc(s) := E[e _sTfe ]. In the following we show that the 
Laplace transform is the only information we need to compute the moments. 
When the sojourn time in state k € E expires, the environment jumps to 
state j € E with probability pkj. Denoting by P := {pkj}k,jeE the rout- 
ing matrix that we assume irreducible and with no loss of generality with 
Pkk = 0, we can define the reverse routing matrix 

(2i) Q:=n -iptn, 

where denotes the transpose of the matrix P, II := diag(7r) a nd tt is the 
statio nary distribution of the Markov chain generated by P (see iBremaud , 



19991 . §6.1). 

We assume that when the environment is in state k € E customers arrive 
according to a Poisson rate Xk > 0. Each of them brings an independent 
request of service, a, that is exponentially distributed with rate [i. All 
servers work at constant speed f3f- = fJ,k/[i> < 1- To avoid trivial cases we 
assume that fi, /?, A > where A := max/j g £ Xk and (3 := max^g (3k- 
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By the results in iD'Aurial (120071 ) the system is stable and we are allowed 
to study its stationary regime. 

We then look at the system at time and we count the number of 
customers still in the system. We order them according to their arrival 
times {uh}hez with Uh < u^+i and u_i < < no, and we denote by 
G(a) := 1 — e _/WT , a > 0, the common exponential distribution function of 
the {(Jh}h&- 

The /i-th customer, h < 0, will be in the system at time iff its service 
time, ah, is bigger than the work done by the server it has occupied dur- 
ing the time interval [uh,0). We denote this quantity by i*r(u/j) and, as 
the subscript shows, it is a random quantity that depends on the random 
environment T. Its value can be computed in the following way, 

(2.2) F r (u) := f f3 m dt, u < 0. 

J u 

Denoting by N the number of customers in stationary regime we have 
that it is given by 

(2.3) N = J2 1 Wh>F T (u h )}, 

h<0 

where 1{-} is the indicator function of the set {•}. It is helpful to rewrite the 
numerable collection of indicator functions appearing in expression (|2.3p in 
the following equivalent way 

IK > F T (u h )} = 6 {Uhiah) (A v ) 

where <5( U)CT ) denotes a Dirac delta measure with center (u, <r)elx M + and 
the set A-p C M~ x IR + is given by 

(2.4) A r := {(u, a) € IR~ x M. + : a > F r (u)}. 

This alternative formulation allows the decoupling of the sequence {(uh, ah)}hez 
and the function Fr(u) both depending on the realization of the environ- 
ment r in the computation of the quantity N. Indeed we can express the 
stationary number of customers in the system in the following way 

(2-5) ^ = EW)( A r) = W, 

h<0 

where 

(2-6) ^r:=£<W) 

is a point process which locates one Dirac delta measure at each arrival 
point {(uh, <T f,.)\h(=7.- For the theoretical b ackgrou nd and definiti on of point 
processes see 



Dalev and Vere-.TonesI (Il988l . §7.1) or lResnickl (|l98l §3.1). The 



subscript T stays to denote that A/"r depends on the random environment 
by the sequence of arrival times {tihjhez- Indeed given a realization 7 of 
the process T, the sequence ({uh}h£z\F = 7) belongs to an inhomogeneous 
Poisson process with intensity rate XjUa, u £ R. By Proposition 3.8 in 
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Resnickl (| 19871 ) it follows that J\f r \T = 7 is still a Poisson process, now on 
1R x R + , with intensity measure 



X 7 (A) :=E[Af r (A)\T = i\ 



X y ^duG(da), icRx 



Finally A/r is a doubly stochastic Poiss on process or, more briefly, a Cox 



process (see lDalev and Vere-Jonesl . Il988l . §8.5), i.e. a Poisson process whose 



intensity measure is itself random and given by 

(2.7) X r (A) := E[M T (A)\T] = [ X r(u) duG(da), AcRxR + . 

J A 

It is well kno w that the fidi distributions o f a Cox process are of mixed Pois- 
son type (see lDalev and Vere-Jonesl . Il988l . Corollary 8. 5. II), or equivalently 
that for any set A C R x R + 



(21 



Mr(A) ~ Po(|A| r ) 



is a Poisson random variable whose parameter is itself random with value 
\A\r = K[Mr(A)\T]. \A\r can be geometrically interpreted as the measure of 
the set A according to the measure Ar(-), i.e. \A\r = Ar(^4). 
From relations (|2.5p and (|2.8p we finally get that 

(2.9) N ~ Po(|^r|r), 

a mixed Poisson random variable with random parameter |^4r|r = Ar(^4r)- 
Figured] shows an example of realization where the random environment has 



A r 




Figure 1 . Example of realization. 

K = 5 states: the dots are the centers of the Dirac deltas of the point process 
A/r, while the piecewise linear function Fr(u) denotes the lower bound of 
the set of integration Ay ■ The customers present in the system at time are 
then the ones whose dots fall in the set ^4r; in the shown example N = 2. 
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Example 2.1. The easiest case is when the environment process is constant, 
K = 1, so that the system reduces to a classical M/M/oo queue. In this 
case the set Ay is deterministic, given by {(it, a) : u < 0, a > (3\u\)}. From 
equation (|2.7j) we get 

/ AG(At) du = - 1 - G(u) = -E[<r] = — 

-oo 7ct>/3|u| P Jo P PH 

and we obtain the known information that N ~ Po(-^), i.e. the stationary 
number of customer in the system is Poisson distributed. 



3. Computing the factorial moments 

Before beginning to compute the factorial moments of the random variable 
N, it is worthwhile to review some basic results about the different kinds of 
moments and their relations with the various generating functions. A good 
reference a bout the following rela t ions e specially in connection with point 
processes is Daley and Vere-Joned ( 19881 ). Chapter 5. 



Given a random variable X, we denote by ipx(s) '■= lE[e s "^] its moment 
generating function and by 4>x{z) := E[z^] its probability generating func- 
tion. 

The factorial moment of order i of X, /J^ is defined as 



n-pr, 



n=0 



where p n = Pr{X = n} and n- := n(n— 1) • • • (n—i+1) is the falling factorial. 
It can be directly computed by the i-th derivative of the probability generat- 
ing function, i.e. = lim 2 ^i (p^ (z). Knowing the factorial moments of X 
it is then easy to compute its moments, in the sequel called raw moments to 
distinguish them from the factorial ones. Indeed, by taking the expectations 



on bo th sides of the following known equivalence (lAbramowitz and Stegun . 
1964 ) 



i=0 

where is a Stirling Number of the Second Kind, we obtain the following 
relation between the n-th moment of X, := E [X n ] with m5 := 1, and 
the factorial moments of order i < n, 



(3-1) = 

o 



The reverse relation is obtained by using the Stirling Numbers of the First 
Kind, (see Abramowitz and Stegun . 19641 ). that satisfy the following 
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known relation 

x^ = ]T Sj {n) x n , 

n=0 

so that, taking the expectations of both sides, finally we get 
(3.2) /» = £X B >m£>. 

n.=0 

It is interesting to notice that relation (13. lj) comes directly from using the 
fact that tpx{s) = cj)x{e s ) and that = lim s ^o ip^\s) . Indeed, 

m n 

Hm#( S ) = = £6»<$(1), 

i=0 

where in the last equation we used Faa di Bruno's form ula for the expansion 



of de rivatives of order n for composition of functions (see lAbramowitz and Stegun . 
19641 ) and the fact that lim s ^ = 1. 



A random variable X is called mixed Poisson when there exists a 
non-negative random variable Y such that X = Po(Y), or equivalently 

(X\Y = y) — Po(y), where the operator = denotes equality in distribution. 
In the case X were a mixed Poisson random variable we would have that 

4>x{z) = ^ Y {z - 1), 
so that taking the derivatives of order n, we get 

lim <f>x\z) = lim iPy\ z — 1) — hm tfjy\s), 

or, in other words, that the factorial moments of X are directly the raw 
moments of Y, 

and the latter often are easier to compute. 

This is exactly what happens in our case where, as shown by relation (|2.9p . 

is a mixed Poisson and that is why we are interested into its factorial 
moments rather then directly its raw moments. Indeed we have that the 
following relation holds 

(3-3) fN = m \Z\r' 

and our task reduces to the computation of the raw moments of the measure 
of the random set Ay . 

4. Computing the raw moments of |^r|r 

In this section we compute the raw moments of the measure of the set 
A-p, defined in (12. 4j) . when measured by the random intensity measure Ar, 
defined in (|2.7j) . We use a fixed point technique and to this aim we look at 
a modified environment process, Tq, that is the Palm version of the process 
T, i.e. we assume that at time it has a transition. We denote by k € E 
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the last state it has assumed before 0, i.e. k := ro(0~), and by T k its 
corresponding sojourn time. While, as depicted in Figured! for the process 
r the sojourn time in the last state before would be given by a residual 
sojourn time, for the process To it is distributed as any other sojourn time 
corresponding to the same state. We define by A^ := (^4r |ro(0 — ) = k), 
k £ E, the set Ar conditioned to the event that the last state occupied by 
the environment before is the state k, and we call |^4ofc| its measure, i.e. 
\A ok \ := (Ar o (Ar o )|r o (0-) = A;). 



V)3 



, 3 i,Ao5 



T 3 I T, I T 5 



'03, 




&T 2 

"/34T4 
_/3 3 T 3 

/3iT, 

> 5 T 5 
P3T3 



Figure 2. Decomposition of Aqz as C03 U Tr 3l p 3 T 3 Ao5. 



Figure [2] shows an example of the set Ao k when k = 3, together with 
its decomposition in the set Cofc and the set 7y 3 ^ 3 t 3 Aq^. To this we have 
defined by Cofc the restriction of the set Aqi~ up the last transition of the 
process To before time 0, i.e. 

(4.1) C ok := A ok n {(x,y) € R~ x R+\\x\ < T k }, 
\C k\-= (Aro(Co fc )|r o (0-) = fc) 

and by T s ,tA the (— s, t)-translated version of the set A, i.e. 

(4.2) Tstt A := {(x,y) &R 2 \(x + s,y - t) € A}. 

We denote by j £ E the state of the environment before the last transition 
before time 0, i.e. j := To(— T^T), so that, — T k being a regeneration point 
for the process To, we have the independence of the sets Co k and 7r fe ,/3 fc T fc A)? 
conditioned to the values of the states before and after the transitions, i.e. j 
and k. f3 k T k = Fr (—T k ) is the exact amount of work the non-empty servers 
have done during the time interval [—T k , 0) being in state k. 

By noticing that the set 7r ki /3 k T k Aoj has measure equal in distribution 
to \T 0ATk A 0j \ := {Xr (%,i3kT k Aoj)\T (0~) = j), we can write down the 
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following set of stochastic equations 

K 

(4.3) \A ok \ 4 |C ofc | +Y, 1 i k ^ j}\ T 0AT k A 0j \, 

i=i 

where the indicator function \{k <— j} selects the backward state transition 
of the environment from the state k to the state j; this would happen, 
according to definition (|2.ip . with probability q k j. 

Thanks to the fact that along the vertical axis the measure function is 
given by G that is exponential we have that the following result holds: 



Lemma 4.1. Given the transformation T Si t, defined in |^TI|), we have that 

(4.4) |T ,^ r |r = e" M Vr|r, 

for any random set A-p : T — ► B(M~ x R + ). 
Proof. By using definition f|2.7j) we have 

A r (T 0it A r )= / Xr^due-^da 
Jr , t A r 

= e-^ [ A r(u) du e-^ a - l) da = e~^ \p (Ap ) . 



□ 



By using Lemma 14.14 equation (|4.3p simplifies in the following 

K 

(4.5) \A ok \ = \C ok \ + <T^ kTh ^ k <- filAvl* 

j"=i 

that is the starting point to prove the following main result: 

Theorem 4.2. Let us define fn^ € M. K as the column vector whose k-th co- 
ordinate is the n-th moment of the random variable \Ao k \, i.e. m^' := 
then the following relation holds 

(4.6) ^(- 1 ) l f n ) Rn " lB -^ )=0 ' 
i=o w 

where R := diag (p k ), Pk '■= X k /p k o-nd the matrix~B n : = diag (r^ 1 (np k )) — Q. 
The matrix B n , n > 0, is invertible and therefore it is possible to express the 
n-th moment vector friQ in terms of the previous ones, w!q , i = 0, . . . ,n—l, 
in the following way 



n-l , v 

(4.7) mf' = E(- 1 ) n " 1 " i ( n )B^R-B n m« 
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Proof. We first compute the values of the variable | Cofc | in the following way 

r-T fc 



IC, 



Ok 



dx = p fc (l - e"^ Tfc ). 



Then substituting its value in equation (|4.5p . it gives 

A 

(4.8) |A ofe | = p fc (l - e^ Tfc ) + e"^ ^ l{k <- j}\A 0j \, 

i=i 

that can be rewritten as 

A 

(4.9) \A ok \ - Pk ^Y, l i k ^ j'KIAyl " Pfc)e- MfeTfe - 

i=i 

We denote by ipok(s) '■= E [e s ' j4 ° fc '] the moment generating function of |^4ofc| 
so that applying the exponential function to both members of equation (|4.9p 
previously multiplied by s and then taking the expectation, we obtain 

A 



' A 



Last expression can be written in matrix form in the following way 



(4.10) 



sR 



Ms) = E 



sRT 



where T := diag(e _/XfeTfe ) and where with notation v(W), with W a diagonal 
matrix, we denote a vector whose k-th component is v^Wkk)- We use then 
the following matrix formulas for derivatives 

(4.11) D^[e" sW v(s)] = ^(-l) n " i (' nN )e" sW W n ^Z)«p(s)], 
and 

(4.12) D^[v(sW)]=W n v^{sW), 

to compute the ra-th derivative of both sides of equation (|4.1U|) so that 

^(-l)-^ n V-^R--#)( s ) = 
i=0 W 



E 



E 



i=0 



i=0 



-R"-i T »(Q^W)( sT ) 
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Remembering that mg"' = lim s ^o ( s ) an< ^ taking the limit of last ex- 
pression as s — > 0, we get 

(4.13) ^(-lf-^^R^m? =E Q-l)"- i MR' 1 -TQ4 i) 



i=0 



Multiplying on the left side by (-l)~ n E[T n ]~ 1 , the last expression can be 
easily rearranged in 

(4.14) f^ir^V-pEfTT 1 " Q]m? = 0, 

that gives the result. The invertibility of the matrix B n for n > comes 
from Lemma I A. 21 □ 



It is remarkable that it is possible to express equation (|4.6p in terms of 
the forward transition chain P. The result is contained in the following 
corollary whose proof comes from simple matrix computations. 

Corollary 4.3. A result similar to equation ( L/.6p is valid for the row vector 
^t(n) ._ /-.(») 



m j , that involves the matrix P instead of the matrix Q, i.e. 

(4.15) f>l)«("WnB' B R»-* = 0, 

i=0 ^ ' 

where := diag (r / 7 1 (?i/i/ c )) — P. TTte matrix B^ zs non-singular when 
n > 0. 

Given the raw moments of the |Ap |, we can successively compute the 
moments of the measures of the sets := (Ar|r(0) = k), k £ E. Following 
previous definitions we define := n^f^j- Similarly to equation (14. 3ft we 
have the following equation 

K 

(4.16) |^| 4 + ^ l{fe - i}|T OAT .A oi |, 

i=i 

with |C?| = Pfc(l — e _AtfeT fc). T£ refers to a residual sojourn time of the 
environment in state k € E. We define by t£(s) the Laplace transform of 
the distribution function of T£ and it is related to the one of T^, Tfc(s), by 
the relation Tjt(s) = ffc(l — Tfc(s))/s, with ffe := E[Tfc] -1 . For the vector of 
raw moments m^ n ^ the following theorem holds. 



Theorem 4.4. The vector rnS n ^ satisfies the following relation with the vec- 
tor fn^ 

(4.17) f^-l^W-^W - E„4 i} ] = 0, 

i=0 
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with E n := diag(r^(n//A;)/T)fc(n/Xjfc)). Therefore the vector can be com- 
puted from the previous moments {tfi"}i< n and the corresponding vectors 
{Aq }i<n in the following way 

n-l , v 

(4.18) # =E n m( n) + J^(-l) n - 1 - i ( n R n ~ l [m« -E n ra[ ], 

i=o W 

finally 

(4.19) mW := = m^vf. 

Proof. Starting by equation (|4.16j) and following the same calculations that 
brought us from equation f|4.3|) to equation (|4.14j) . we get 



(4.20) ^(-lfWR^pStT*"]" 1 ^ - Qm?] = 0, 

t=o 

that after subtracting equation (|4.14j) gives 

(4.21) ^(-l^W^E^r 1 ™® -EpTT 1 ^ ] = 0. 

and by multiplying on the left by E [T* n ] we finally get the result. □ 



In order to check our re sults we compare equati on ( 14. 151) for the expo- 
nential case with results in O'Cinneide and Purdue! (jl986l ) here repeated in 

formula (I4.22p . For this case since I? = we have that = Aq . 



Remark 4.5. It is worth to notice that in lO'Cinneide and Purdue! (| 19861 ). 



they actually computed the factorial moments of the random row vector 
(N 1{Tq = k})k£E while here we compute the factorial moments of the row 
vector (iV|ro = k)k£E- This explains the presence, in formula (|4.22[> . of the 
additional factor given by matrix II. 

Corollary 4.6. In case the sojourn times are exponentially distributed 
with parameters ff. we have that with n > 

(4.22) (mt( n) n) (nM - G) = n (m t(n_1) II) A 

where M := diag(/ifc), A := diag(A^) and G := T(P— I), with T := diag(ffc), 

is the generator of the Markovian Environment. 

Proof. When the sojourn times are exponentially distributed we have 
that 

E [ T nj-i = ^MT" 1 + 1, 

that implies 

TB n = nM - H, 

with H := T(Q — I) being the generator of the reverse-time Markov process. 
By multiplying on the left both sides of equation (14. 6ft by T and noticing 
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that it can commute with the powers of the matrix R, we can rewrite (|4.6j) 

as 

(4.23) ^(-l) n " i ( n \ R n ~>M - H] m« = 0. 

By Lemma EI with U = M, D = R, V„ = H and u (n) = mW, and 
imposing = 1, we get that the unique solution of (|4,23p is given by 

(nM - H)m {n) = nRMm^" 1 ), 

that transposed reduces to 

^(nr^nM - H*) = nm^^A. 

Multiplying on the right side by II and simplifying we get 

fflt("'n(nM - n^II) = nm^ n - l) UA, 

that gives the result after noticing that G = II _1 H^II. □ 

5. Some explicit formulas - Case K = 2 

Formulas (|4.7p and (|4.18p show that generally to find the n-th moment 
of the random number of users in the system involves in a complex way the 
knowledge of all prev ious moments. Reversel y the exponential case, that 
was already solved in lO'Cinneide and Purdue! (|l986l h is easier as the n-th 



vector of moments is related only by a factor to the (n — l)-th one. That 
was anyway hidden in a non-trivial way in formula (|4.7|) so that there could 
be some other special cases where an easier expression holds. 

In this section we have a look to the case when the environment has only 
two stages, i.e. K = 2. 

This is a very special case and when the sojourn times are all assumed 
exponentially distributed, it is known how to compute the complete dis- 
tribution of the numbe r of customers in the system at ste ady state (see 
Keilson and Servil (|l993MBavkal-Gursov and Xiao! (120041 ) and lD'Aurial l|2007l )). 



We give for this case an explicit formula to calculate the factorial moments 
in terms of the Laplace transform of the sojourn time in state 1, when the 
sojourn time in state 2 is exponentially distributed. 

By rewriting in more explicit form equation (|4.9p for the case K = 2 we 
get 

(5.1) L4oi I- Pi = (14)2 1 -pi)e- mTl 

(5.2) |4 2 |-p 2 = (!4i!-p 2 )e" At2T2 . 

We define := E[(|4fc| — Pi) n ] an d take the mean of the n-powers of 
expression (|5 ip so getting 

(5.3) m ( $=m ( $T 1 {n» 1 ). 
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By adding and subtracting p\ to both sides of equation (15. 2p we get 

(|Ab| - Pi) - P* = (I An I - pi)e~^ - p*e~^ T \ 

with p* = pi — p\. Then using equation (|5.1j) we obtain a recursive equation 
involving only |A) 2 | — Pi> 

(|A) 2 | - Pi) - P* = - Pi)e-^ Tl e-^ T2 - p*e-^ T *. 

Taking the n-th power and then the expectation of both sides we get 

Et-ir^Vr^of =r 2 (n M2 )E[((|A 02 | - Pl )e-^ -p*) n ] 

=T 2 (^ 2 )^(-if- i f n Vrvi(^i)4t 

i=0 w 

that, taking into account equation (|5.3p . can be rearranged to get the fol- 
lowing 

(5.4) ^(-lr-f^pr^- 1 ^)^ 1 ^) - l)mg = 0. 

i=0 W 

Theorem 5.1. Assuming that the sojourn times of state 2 are exponentially 
distributed, i.e. T 2 ~ Expifi), the solution of formula |5.^[ j is given by 

(5.5) 4H^)"n ^fe^r i - 

V f 2 / f = i Tl (^l)T 2 (^2) - 1 

and therefore 

'p 2 p*V _i. ,17 « Tf 1 ((i - l)pi) 



(\ n 
T 2 / 



1 7j (ipi)r 2 (ip 2 ) - 1 



(5.7) 4?=E(-)pr^t *=1,2. 

i=0 W 

Proof. Substituting r^ 1 ^) = 1 + s/f 2 in equation (|5.4p and rearranging it, 
we get 

n / \ / \ ~ (*) 

(5.8) E(- 1 ) n "\l)^ ( n T 2 ~ {Tlim) ~ 1} ) ^) = °" 

By applying Lemma [A. II in the scalar case, with U = (p 2 /f 2 ), D = (p*), 

V n = (ri(n/xi) — 1) and = (rh^ /n(npi)), we notice that a set of 
solutions is given by 

~ (n) ~ (n-1) 

(5.9) W) - n(n Ml )) = np/ 2 ^ 



Ti(nm) t 2 n((n - l)pi) ' 
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that is then uniquely defined given that m^' = 1. Therefore equation (|5.5p 
holds. Equation (|5.6p results by applying (|5.3p to (|5.5p and finally equation 
(I5T71) comes from the fact that m$ = E[((\A ok \ - pi) + p 1 ) n ]. □ 

Example 5.2. Case T\ ~ Exp(-fi). In this special case equation (|5.5p sim- 
plifies in 

~{n) _ n (n/Mlf 

01 " P * (ri/Mi + r 2 //, 2 + 1)"' 

with z n := • • • (i+n—1) being the rising factorial (in lAbramowitz and Stegun 



(|l964l ) it is denoted by (i) n ). Therefore the moment generating function of 



Am | — pi is given by the Kummer fu nction M(ti /pi , Ti/ fjii+rz/ fJa + 1, p*s) 
(see lAbramowitz and Stegunl (19641 )). i n accord a nce t o what is shown in 



Bavkal-Gursov and Xiad 42004J) and in iD'Aurial (|2007l l [ in there, it is de 



noted by ON (-«)]■ 

The following example is a new result that generalizes the one of Example 
1531 



Example 5.3. Case Ti ~ Gamma(K,Tj x ). For this case we have that 
r^ 1 (kpi) = (1 + kp\jf\) K . Therefore equation (|5.5|) simplifies in 

-in) = n!p?[(n/ m )T 

01 niLi Kn/m + iY(j2lv2 + - (n/w)^//^)] ■ 
6. Conclusions 

In this paper we showed that using a matrix-geometric approach it is 
possible to solve the problem to find the factorial moments of the random 
number of customers in an M/M/oo system when its parameters are mod- 
ulated by a semi-Mar kovian random environment. We showed that this is 
possible by looking at this random variable as the random measure of a bidi- 
mensional random set by a mixed Poisson process. Finally the case when the 
environment has only 2 states is more deeply investigated and it is shown 
that explicit formulas are obtainable given that one state has exponential 
sojourn times. It is then plausible to believe that for this last case it would 
be possible to get an explicit expression for the complete characteristic func- 
tion. 

Appendix A. Technical Lemmas 

Lemma A.l. Given the matrices U, D, {V n } n >o such that for any < n < I 

the matrix (nU — V n ) is invertible than the system of equations 

(A.l) ^(-lf-^W-^nU- Vi]v® =0, 0<n<l 

has a family of solutions v ( n \ < n < I, given by 
(A.2) (nU - V n )v^ = n DU v {n ~^ , 
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that reduces to a unique solution once given the vector that has to satisfy 
the relation Vqv = 0. 

Proof. The proof of the lemma is immediate once we prove that a set of 
vectors obeying to relations (|A.2|) for n < I with Vot^ ) = satisfies as well 
the following equation for < n < I 
(A.3) 

f(n) 



2(-l)'- w QD w [iU- Vi]v® = (-l) l - n (^jB l - n [nU-V r 
We prove it by induction. Assuming n = 1, we have that 



(-l^D^U- V ]£ (0) = (-l) f - 1 D i [ZU] 



(0) 



= (-lj'^ID'-^Ui/M 

= (-i)'- 1 Q)D'- 1 (u-v 1 )i;w > 

so that the base of the induction holds. Now assuming (|A.3|) valid for n < Z 
we have for n + 1 that 

^(-ly-^f^D'-^ZU-VjtfW 

i=0 

= (-l)'- n ^D'- n [nV - Vnjv^ + (-IJ'-^Qd^IIU - V, 

= (— l)'~ 1 ~ n r J D^ n (Z - n)LHT (n) 

= (-lj'-^QD^lI-njDWW 

and by equation (|A.2I) . 

= ^ l^-D i - n - 1 ((n + 1)U - V n+1 )^(™ +1 ) 

\n/ n + 1 

= ( | V^+i)^ + 1)U - V n+1 )<;( n+1 ) 

□ 



The proof of next result follows closely the one of Lemma B.l in ISeneta 
([1983). 

Lemma A. 2. The matrix B n , n > in equation is invertible. 

Proof. By noticing that all the diagonal entries of the matrix diag (t^T 1 (n/j,k)) 
are strictly positive, in order to prove the invertibility of the matrix B n we 
are left with proving the non-singularity of the matrix (I — D n ), with 

D n := diag (r fc (n// fc )) Q. 
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We assume with no loss of generality that the states of the random environ- 
ment are ordered according to the increasing values of so that the states 
with null betas have the lowest indexes. Being (3 > we know that they 
are in number Kq < K. If we compute the I power of the matrix D n we get 
that in the limit it converges elementwise to 

(A.4) T> l n — > ( X fo J V asZ-oo, 



t 

where ttq is the vector containing only the first Kq coordinates of the vector 
7T and 1 is defined as a vector with all coordinates equal to 1 and whose 
dimension depends on the context. 

The above result comes from the fact that Q' — > 1 tt^ as I — > oo and from 
knowing that Tfc(n/Xfc) = 1 when k < Kq and r[(n/Xfc) — > when Kq < k < K 
as I — > oo. For any I > the following relation is valid 

(I - D„)(I + D n + • • • + D^ 1 ) = I - Di, 

and the matrix on the right side converges to 

1 - D n -» ( I J ^ °°' 

The determinant of the limit matrix in the right side of last relation is equal 
to the determinant of the matrix (I — Tttq) that is positive by applying 
Lemma B.l in Seneta ( 198ll ) to the strictly s ubstochastic m atrix 1 ttq. 



Following the reasoning in Lemma B.l in ISenetal (|l98ll l. as the deter- 



minant is a continuous function in the space of matrices with elementwise 
convergence we get that there exists some I > such that the determinant 
of the matrix (I — DjJ is positive. Therefore it follows that the product of 
the determinants of the two matrix factors in the left side of equation (|A.4|) 
has to be positive, that concludes the proof. □ 
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